#ifndef WARRAY_IMPL_WARRAY_DFT_H_INCLUDED
#define WARRAY_IMPL_WARRAY_DFT_H_INCLUDED

#include <complex>
#include <cmath>
#include "../warray.h"

namespace warp {
namespace detail {
namespace {
const double pi=3.14159265358979323846;
} // namespace

//! \brief Naive implementation of a 1D DFT
template<typename T>
class naive_1d_dft_implementor
{
public:
	//! \brief Initialise a DFT for n elements
	naive_1d_dft_implementor(size_t n) :
		m_coeffs(warp::extents[n][n]),
		m_buffer(warp::extents[n])
	{
		warp::arraynd<T,2> theta(m_coeffs.dims());
		for(unsigned int i=0;i<n;++i)
			for(unsigned int j=i;j<n;++j)
				theta[i][j] = theta[j][i] = T(i*j);
				
		T factor = T(-2.0*pi/n);
		theta *= factor;
		warp::cis(theta,m_coeffs);
	}
	
	//! \brief Compute the 1D dft *not-in-place*
	void fwd(const std::complex<T>* l, int lstride, std::complex<T>* r, int rstride)
	{
		warp::blas::gemv(
			m_coeffs.dim(warp::dim1),
			m_coeffs.dim(warp::dim2),
			m_coeffs.c_array(),
			m_coeffs.stride(warp::dim1),
			l,lstride,
			r,rstride);
	}
	
	//! \brief Compute the 1D dft *in-place*
	void fwd(std::complex<T>* r, int rstride)
	{
		fwd(r,rstride,m_buffer.c_array(),1);
		warp::blas::copy(dim(),m_buffer.c_array(),1,r,rstride);
	}
	
	//! \brief Determine the size of DFT computed
	size_t dim() const
	{
		return m_coeffs.dim(0);
	}
	
protected:
	warp::arraynd<std::complex<T>,2> m_coeffs;
	warp::arraynd<std::complex<T>,1> m_buffer;
};



//! \brief Naive implementation of an N-dimensional DFT
template<typename T, size_t N>
class naive_nd_dft_implementor
{
public:
	//! \brief 
	naive_nd_dft_implementor(const warp::array_extents<N>& dims) :
		m_n_minus_one(dims.slice_one()),
		m_implementor(dims.dim(0))
	{
	}
	
	//! \brief Compute an N dimensional DFT
	template<typename Tptr>
	void fwd(const const_arraynd<std::complex<T>,N,Tptr>& in, const arraynd<std::complex<T>,N>& out)
	{
		bounds_check(in,out);
#ifdef BOUNDS_CHECK
		assert(in.dim(warp::dim1) == m_implementor.dim());
#endif
		for(size_t i=0;i<m_implementor.dim();++i)
		{
			m_n_minus_one.fwd(
				in.slice(warp::indices[i]),
				out.slice(warp::indices[i]) );
		}
		
		size_t M = in.size()/in.dim(warp::dim1);
		std::complex<T>* optr = out.c_array();
		std::complex<T>* oend = optr+M;
		for(;optr!=oend;++optr)
			m_implementor.fwd(optr,out.stride(warp::dim1));
	}

protected:
	naive_nd_dft_implementor<T,N-1> m_n_minus_one;
	naive_1d_dft_implementor<T> m_implementor;
};

// Specialisation for the one dimensional case
template<typename T>
class naive_nd_dft_implementor<T,1>
{
public:
	naive_nd_dft_implementor(const warp::array_extents<1>& dims) :
		m_implementor(dims.dim(0))
	{
	}
	
	template<typename Tptr>
	void fwd(const const_arraynd<std::complex<T>,1,Tptr>& in, const arraynd<std::complex<T>,1>& out)
	{
		bounds_check(in,out);
#ifdef BOUNDS_CHECK
		assert(in.dim(warp::dim1) == m_implementor.dim());
#endif
		// Check for memory aliasing
		if(in.c_array() != out.c_array())
		{
			m_implementor.fwd(
				in.c_array(),in.stride(warp::dim1),
				out.c_array(),out.stride(warp::dim1));
		}
		else
		{
			m_implementor.fwd(
				out.c_array(),out.stride(warp::dim1));
		}
	}
	
protected:
	naive_1d_dft_implementor<T> m_implementor;
};

} // namespace detail
} // namespace warp

#endif

